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Diffusion of self-propelled particles in the presence of randomly distributed obstacles in three 
dimensions is studied using molecular dynamics simnlations. It is fonnd that depending on the 
magnitude of the propelling force and the particle aspect ratio, the diffusion coefficient can be a 
monotonically decreasing or a non-monotonic concave function of the obstructed volume fraction. 
Counterintnitive enhancement of the particle diffusivity with increasing the obstacles crowd is shown 
to be a combinatory effect of the self-propelling force and the anisotropy in the shape of the particle. 

Regions corresponding to monotonic and non-monotonic dependence of the particle diffusivity on 
the obstacle density in propelling force-aspect ratio plane are specified theoretically and nsing the 
simulation results. 
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Objects under simultaneous influence of a self- 
propelling force and a randomly varying force exerted by 
environment have been of great interest in recent years. 
Microscopic organisms m , artificial microswimmers n- 
synthetic motile objects from molecular scale to mi¬ 
crons surface-active colloidal particles [mils] and 
millimeter-sized manmade particles [14H22j are examples 
of such objects. Combination of the two mentioned forces 
has been shown to cause a variety of regimes in dynamics 
of such objects [131 US]. 

The shape of a particle dispersed in a fluid strongly 
affects its diffusion dynamics. An anisotropic particle 
experiences different values of friction coefficient when it 
moves in different directions and the coupling of transla¬ 
tional and rotational motions makes understanding and 
visualization of its diffusion noticeably difficult [331 ■ 
propulsion amplifies the difference between dynamics of 
isotropic and anisotropic particles. 

Obstructed diffusion - diffusion of particles in the pres¬ 
ence of fixed obstacles - is a ubiquitous phenomenon in 
nature. The diffusing object in the presence of obstacles 
could be a simple spherical particle, a rigid anisotropic 
particle or a long flexible polymer. Obstructed diffusion 
of passive particles of all the mentioned types has been 
studied extensively |3SH33]- Presence of fixed obstacles, 
also, has different effects on the diffusion dynamics of 
isotropic and anisotropic particles. In the latter case, the 
obstruction effect depends on the direction of motion; a 
rod-like particle experiences the least (the most) num¬ 
ber of collisions with the obstacles when it moves along 
(perpendicular to) its axis. 

For self-propelled (SP) particles of anisotropic shape 
(such as rod), an important question is how the crowd of 
fixed obstacles affects their dynamics. In this paper, we 
study dynamics of rod-shaped SP particles in the pres¬ 
ence of randomly distributed spherical obstacles using 
the results of molecular dynamics (MD) simulations. We 
find that depending on the magnitude of the propelling 


force and the particle aspect ratio, the diffusion coef¬ 
ficient of the particle is a monotonically decreasing or 
a non-monotonic concave function of obstructed volume 
fraction. This non-monotonic dependence means that 
in a range of the parameters, diffusivity of the parti¬ 
cle increases with increasing the obstructed volume frac¬ 
tion. Such a counterintuitive effect is shown to result 
from combination of the effects of the propelling force 
and anisotropy in the particle shape. In propelling force- 
aspect ratio plane, regions corresponding to monotonic 
and non-monotonic dependence of the diffusion coeffi¬ 
cient on the obstacle density are specified both using the 
simulation results and theoretically. Self-propelling force 
is found to cause the particle dynamics in the presence 
of obstacles to consist reentrant diffusive regimes in ad¬ 
dition to ballistic and super-diffusive ones. 

Let us consider a rod-shaped particle of length Lrodcr, 
mass M and diameter d = la (with aspect ratio Lrod) 
in three-dimensional space, where a has dimension of 
length. Configuration of the particle in space can be 
specified by position of its center, r = {x, y, z), and three 
Euler angles, </), 9 and if [SS]. Considering drag and ex¬ 
ternal forces (torques) in addition to the random force 
(torque) of the thermal noise, the Langevin equations of 
translational and rotational motions can be written as 
Mvi = -PFi-f- r]}{t) and hui = -yfw* + fVi -f ry[(t) 

where i = Ti,T 2 ,|| (see Fig. and ”t” and ”r” stand 
for translational and rotational, respectively. For a self- 
propelling particle, it is assumed that F) contains the self- 
propelling force as well. In the above equations, Vi = v.Ci 
and LOi = Co.Bi where v and Co are linear and angular ve¬ 
locities of the rod in the lab frame. The noise terms 
follow (r?f(t)) = 0, {il^{t)r]]{t')) = 5kAj5{t - t') and 

Sff = 2'y^kBT, where fc = t, r. For motion of a rod¬ 
shaped particle in a fluid, it is known that 7 ^ = 27 |i 

in which 7 ^ ( 7 y) is the friction coefficient the rod ex¬ 
periences when it moves perpendicular (parallel) to its 
axis |36j . The relation between lab- and body-frame 
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FIG. 1: (Color online) Schematic of a rod-shaped particle and spherical obstacles in three dimensions. Shish-kebab model of 
the rod is used for calculation of its excluded volume interaction with the obstacles (middle). Sample trajectories of spherical 
{Lrod = 1) and rod-shaped {Lrod = 11) particles (left and right). The most diffusivity of the rod corresponds to its obstructed 
SP diffusion. 


views of any vector A reads A^ody = 5ft(((), 0,'0 )A;q;,, in 
which 5ft is the rotation matrix [^. Regarding that 
5ft has a singularity at ft = 0, Cayley-Klein parameters 
can be used; qi = sin(|) cos(^^), q 2 = sin(|) sin(^^), 
qz = cos(|)sin(^^) and q^ = cos(|) cos(^^). These 
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parameters follow equations q^ = ^ and ^ = 0 

i=l i=l 

[351I3ZH3S]- Direction of the particle in space can be up¬ 
dated as qi{t+At) = qi{t)+qi{t)At, f = 1, 2, 3,4 in which, 
ft = and Q{q) is a matrix formed by g^ as: 
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To study obstructed diffusion of the particle, the spher¬ 
ical obstacles were modeled by fixed non-overlapping 
spheres of diameter d = la, randomly distributed in¬ 
side the simulation box. Excluded volume interactions 
between the rod and the obstacles were modeled by 
shifted and truncated sphere-sphere Lennard-Jones po¬ 
tential (see Fig. [^. The simulation box was a cube of 
edge length L = 20a and periodic boundary conditions 
were applied. Second order Runge-Kutta algorithm was 
used for integration of the equations of motion [?D] . The 
temperature was kept fixed at T = As the fluid was 
modeled implicitly by stochastic and damping terms in 
the Langevin equation, hydrodynamic interactions were 
not considered. MD time step was r = O.OOlrg, in which 

tq = ^is the MD time scale. The values of the 
friction coefficients were taken as 7 y = Lrod (in units of 

, 7i = 27 ^ and 7 I = M- For given 

values of Lrod, magnitude of the self-propelling force (F^p 
in units of F^) and the obstructed volume fraction ($), 
the simulations were repeated with 10 different realiza¬ 


tions of the obstacles and the averages were calculated. 
$ is defined as the fraction of the simulation box forbid¬ 
den for a spherical particle of diameter a due to excluded 
volume interactions. 

By simply looking at the trajectories of particles of 
different Lrod at different values of Fgp and $, an inter¬ 
esting phenomenon can be seen. With some values of Fgp 
and Lrod and in a given duration of time, the SP particle 
covers considerably longer trajectory in the presence of 
obstacles relative to the space without obstacles (see Fig. 

right panel). This is a counterintuitive phenomenon 
as it is expected the obstruction to decrease the particle 
diffusivity. Sample trajectories of a sphere of diameter 
d = la {Lrod = 1) and a rod of Lrod = H obtained 
from simulations of time length 3 x IO^tq are shown in 
Fig. [3 As it can be seen, diffusivity enhancement due to 
the presence of obstacles doesn’t happen for SP spherical 
particle showing that anisotropy in the particle shape is 



FIG. 2: (Color online) Reduced diffusion coefficient of parti¬ 
cles of aspect ratio Lrod versus il? for different values of Fsp. 
Relative error on the data points is ~ 10 percent. 
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a key parameter in obstruction-induced enhancement of 
the diffusivity. 

As the main measure of the particle diffusivity, we cal¬ 
culate its diffusion coefhcient from time dependence of 
the mean-squared displacement (MSD) and then aver¬ 
age over obstacles realizations. In Fig. dependence 
of the reduced diffusion coefficient, on $ is shown 
for different values of Fgp and Lrod {Dq is the diffusion 
coefficient in the absence of the obstacles). As it can be 
seen, for large enough magnitude of the propelling force 
and the particle aspect ratio, the diffusion coefficient is 
a non-monotonic concave function of <i>. By addition of 
obstacles to the environment, the particle diffusivity goes 
beyond that in the environment empty of obstacles and 
reaches to a maximum and eventually decreases to zero. 
For the spherical particle however, D is observed to be a 
decreasing function of <&, regardless of the value of Fgp. 
In fact, both anisotropy in the particle shape and self¬ 
propulsion of enough strength are needed the enhanced 
diffusivity to occur. 

In addition to the translational diffusion coefhcient, 
we calculated orientational relaxation time of rod-shaped 
particles, t^, by looking at their orientational diffusion 
and from time dependence of (ey(t).e||( 0 )) in the lab 
frame. It is found that the particle orientation relaxes 
exponentially as exp(—t/rj.) from which Tr can be cal¬ 
culated. Also, such orientational relaxation means that 
{Fsp{t).Fsp{0)) = Fsp^(e||(t).e||( 0 )) = exp (-t/r^). 

Considering that the values of Lrod and Fgp deter¬ 
mine the reduced diffusion coefhcient to be monotonic 
or non-monotonic function of $, a question is how re¬ 
gions of Fgp — Lrod plane correspond to the two kind 
of dependencies. The simulation results for answering 


this question are summarized in Fig. ^ To describe 
the simulation results shown in Fig. ISTwe hrst calcu¬ 
late the diffusion coefhcient of a SP rod in the absence 
of obstacles, analytically. To this end, the Langevin 
equation for the translational degrees of freedom at long 
enough time scales where the inertial term is washed out, 
should be solved. In this regime, Vi{t) = -I- 

li 

Fsp{t)Si^\\) where i =||,_Li,_L 2 - Therefore, Axi{t) = 
Xi{t) — Xi{0) = and {Axi{t)Axj{t)) = 



lo lo + {F,p{L)F^p{F))d,^llSj^l\]. 


Considering {Fsp(t')Fsp{t")) = Fsp'^ exp{-\t' - t"\/Tr) 
and MSD = X]i=|| _Li ±2 = 6Dt, D can be ob¬ 
tained as D = _|_^v From the above equation 

-V ^ 


37n 


and that Tr = 


37,) 

= 3'll 
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2kBT 


'^11 F ^ 
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( 1 ) 


The hrst term of D in Eq. which is a decreasing 
function of Lrod, corresponds to the passive diffusion of 
the particle. Its dependence on 7 ^ comes from relations 

D = = 27 || [M]. The second term of D 

in Eq. [^is the contribution of self-propulsion. Sum of the 
two terms shows that F is a non-monotonic convex func¬ 
tion of 7 y = Lrod with a minimum at L*^^ = 2 -^/ 6 ^^ 
as shown in Fig. In this hgure, the simulation results 
for $ = 0 (for comparison with Eq. and three nonzero 
values of $ are shown. Also, dependence of D on F^p at 
$ = 0 and two nonzero obstacle densities are shown in 
Fig. a Non-monotonic dependence of active particle dif- 



Lrod 

FIG. 3: (Color online) Regions in Fsp—Lrod plane in which the 
diffusion coefficient of the SP rod is a monotonically decreas¬ 
ing or a non-monotonic concave function of obstrncted volnme 
fraction. The solid red line shows the theoretical prediction of 
the boundary between the two regions; Fsp = 2-^/6 and 
the circle symbols (connected by dashed line as a guide for 
eyes) are the simulation results for the same boundary. 



Lrod 

FIG. 4: (Color online) Diffusion coefficient of a SP particle 
{Fsp = 3) versus its aspect ratio. Solid line corresponds to Eq. 
a and the simulation results are shown by symbols. Insets: D 
versus Lrod from Eq. unlike the passive diffusion for which 
D is a decreasing function of Lrod, with nonzero values of Fsp, 
D is a convex function of Lrod (up). Dependence of D on Fsp 
for a rod of Lrod = 11 (down). 
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fusivity on its size has also been reported and discussed 
in Ref. [13]. 

The boundary between the two regions of the Fsp—Lrod 
plane shown in Fig. can be described by considering 
the non-monotonic dependence of D on shown in 

Fig. a It seems that this boundary corresponds to the 
relation between Fgp and L^od that minimizes D\ Fgp = 
For each given value of Fsp, D has different 
dependencies on the particle aspect ratio for Lrod < 
and Lrod > L*^^ (see Figs. aand^). For values of aspect 
ratio, Lrod < the diffusion coefficient at $ = 0 

is a decreasing function of 7 y. In this region, addition 
of obstacles to the system causes the friction coefficient 
to increase effectively and hence the diffusion coefficient 
to be a decreasing function of <&. For values of aspect 
ratio, Lrod > F*odJ however, the leading term of D is the 
second term in Eq. [l] which comes from self-propulsion. 
This term is an increasing function of 7 y and addition of 
obstacles has a completely different effect on D relative 
to the case Lrod < i*od- For Lrod > Frodi addition of 
obstacles (increasing $ from zero), first causes Tr and 
hence D to increase. Addition of more obstacles to the 
environment eventually causes suppression of the particle 
motion and decrease of the diffusion coefficient. As it 
can be seen in Fig. the theoretical prediction of the 
boundary between the two regions agrees well with the 
simulation result. 

For a sphere and a rod-shaped particle of aspect ra¬ 
tio Lrod = 11, MSD versus time is shown in Fig. 
Passive dynamics of the rod consists of ballistic, sub- 
diffusive and diffusive regimes at short, intermediate and 
long times, respectively. Self-propelling force causes a 
supper-diffusive regime to appear in the particle dynam¬ 
ics at intermediate time scales. In addition, reentrant 
diffusive regimes in the particle dynamics can be seen. 
Similar behavior has been reported for SP particles in 
two dimensions [23) . The diffusive regime at intermedi¬ 
ate time scales is related to the motion of the rod before 
relaxation of its orientation. Long-time diffusive regime 
from which we calculate the diffusion coefficient, corre¬ 
sponds to time scales that the rod forgets its orientation 
and experiences random walk. In Fig. for a rod of 
aspect ratio Lrod = H, TV versus 4) is shown at four dif¬ 
ferent values of Fgp. As it can be seen, is a rapidly 
increasing function of $, Fgp and Lrod- 

In a stochastic model for nucleosome sliding along 
DNA under an external force, non-monotonic dependence 
of the diffusion coefficient on ligand concentration has 
been reported in Ref. [22] . It has been shown in this work 
that non-monotonic dependence of the diffusion coeffi¬ 
cient on the ligand concentration happens only for forced 
diffusions and for large enough magnitude of the driving 
force. In fact, one dimensional nature of the model stud¬ 
ied in Ref. [42] provides an inherent spacial anisotropy 
and makes this model very similar to the motion of a SP 



FIG. 5: (Color online) Log-log plot of MSD versus time for 
obstructed diffusion of a passive and a SP rod of Lrod = 11 
and a SP sphere at >I> = 0.4. Inset: Orientational relaxation 
time versus <1? for a rod of Lrod = 11. 

rod obstructed by fixed particles. In another model for 
study of the driven diffusion of particles in the presence 
of obstacles, non-monotonic dependence of the diffusion 
coefficient on the driving force has been observed 1231- 
I2H] . General similarity between the mentioned models 
and the problem we studied here originates from the fact 
that in all of them, two factors namely the driving/self¬ 
propelling force and the obstruction compete with each 
other. 

In summary, we have shown that obstruction may en¬ 
hance the diffusivity of SP particles. Depending on the 
parameters such as the particle aspect-ratio and the mag¬ 
nitude of the self-propelling force, maximum diffusivity 
can be achieved by tuning the obstacle density. Regard¬ 
ing that the fabrication of SP particles of various shapes 
is possible nowadays, the introduced phenomenon could 
be of interest for drug delivery and many other applica¬ 
tions. 
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